www.gusucode.com > 循环自相关函数工具箱源码程序 > matlab代做 修改 程序循环自相关函数工具箱/cyclostationary_toolbox/cyclic_cross_correlation_fast.m
function R=cyclic_cross_correlation_fast(x,y,T,max_tau) % % CYCLIC_CROSS_CORRELATION_FAST % % calculates the cyclic cross correlation between % two signals x,y at frequency alpha=1/T using a fast % approximation based on the synchronous average of the % time varying autocorrelation. Fundamental signal % period can be defined as a single period or as a sequence % of once per period pulse times. % % R(k*alpha,tau)=E{x(t-tau/2)y(t+tau/2)exp(-jk(alpha)t)} % for k=0 ... 1/alpha % % R=cyclic_cross_correlation_fast(x,y,T,max_tau) % % calculate cross correlation up to max_tau time lags % % if T is a scalar, then T defined the fundamental % cyclic period % % if T is a vector, then T defines a series of once % per revolution impulses % File: cyclic_cross_correlation_fast.m % Last Revised: 23/4/98 % Created: 24/11/97 % Author: Andrew C. McCormick % (C) University of Strathclyde % Simple error checks if nargin~=4 error('Incorrect number of arguments for function cyclic_cross_correlation_fast'); end if T(1)<1 error('Synchronous period must be larger than 1 in function cyclic_cross_correlation_fast'); end % Ensure that vectors are the right sizes and shapes if length(x)>length(y) x=x(1:length(y)); end [rows,cols]=size(x); if rows>cols x=x'; end [rows,cols]=size(y); if rows>cols y=y'; end % Remove excess samples due to non-integer sampling if length(T)==1 % remove jitter samples if non-integer T if T~=floor(T) cp=1;np=1; while cp+T<length(x) cp=cp+floor(T); np=np+T; if (np-cp)>1 x=[x(1:cp-1) x(cp+1:length(x))]; y=[y(1:cp-1) y(cp+1:length(y))]; np=np-1; end end end n=floor((length(x)-2*max_tau-1)/T); else % extract time series correlated with periodic pulses rot_positions=T; T=floor(median(diff(rot_positions))); nx=[]; ny=[]; n=length(rot_positions)-2; for k=1:n; cp=rot_positions(k); nx=[nx x(cp:cp+T-1)]; ny=[ny y(cp:cp+T-1)]; end nx=[nx x(rot_positions(n+1):rot_positions(n+1)+2*max_tau+1)]; x=nx; ny=[ny y(rot_positions(n+1):rot_positions(n+1)+2*max_tau+1)]; y=ny; end % Compute time varying cross correlation r=zeros(2*max_tau+1,floor(T)); temp=zeros(floor(T),n); t=(1:floor(T)*n); for k=-max_tau:max_tau temp(:)=x(t+max_tau).*y(t+k+max_tau); r(k+1+max_tau,:)=mean(temp'); end % Take DFT of time varying correlation with appropriate phase change % to compensate for time shift R=zeros(2*max_tau+1,floor(T)); for k=-max_tau:max_tau R(k+1+max_tau,:)=exp(-i*pi*((0:floor(T)-1)/T)*k).*fft(r(k+1+max_tau,:))/T; end